# ------------------------------------- #
# Replication code for:
#
# Rathbun, Brian, Christopher Sebastian Parker, and Caleb Pomeroy "Separate but Unequal: Ethnocentrism and Racialization 
# Explain the 'Democratic' Peace in Public Opinion," American Political Science Review.
#
# This script reproduces our reanalyses of Robert Johns and Graeme Davies (2012) "Democratic Peace or Clash of Civilizations? Target States 
# and Support for War in Britain and the United States," Journal of Politics.
# ------------------------------------- #

# --- libraries --- #
library(psych)
library(ggeffects)
library(ggplot2)
library(texreg)
set.seed(1912)

# --- set working directory
setwd("~/Downloads/replication_files/")

# --- import data
johns_davies <- readRDS("data/johns_davies.rds")
nrow(johns_davies) # N = 1035, as mentioned in the paper's appendix 

# --- ethnocentrism measure
# reduce ethnocentrism items to unidimensional factor scores
ethno_data <- cbind((6-johns_davies$MISC_w1_q43), johns_davies$MISC_w1_q44, (6-johns_davies$MISC_w1_q45)) # higher values = higher ethnocentrism
ltm::cronbach.alpha(na.omit(ethno_data)) # as footnoted, cronbach alpha = 0.78
ethno_fit <- fa(ethno_data, nfactors=1, rotate="varimax", scores=TRUE, fm="pa")
ethno_fit$loadings # as footnoted, SS loadings = 1.68
johns_davies$ethno_scale <- ethno_fit$scores[,1] 

# also include a simple additive version of ethnocentrism
johns_davies$ethno_scale_additive <- rowSums(ethno_data)

# --- regression models
johns_davies$SPEC_W1_E1casuals <- factor(johns_davies$SPEC_W1_E1casuals)
johns_davies$democracy <- factor(johns_davies$democracy) # democracy is the main treatment effect of interest

# we focus primarily on the interaction between the democracy treatment and respondent ethnocentrism
# we further consider interactions with authoritarianism and social dominance orientation, each of which might
# also amplify hawkish tendencies
summary(mod_full <- lm(SPEC_W1_E1q1 ~ democracy + Christian + SPEC_W1_E1casuals + socialdom + 
                         Nationalism + authoritarian + gender + ConID + LibID + milpac + age65plus + 
                         ethno_scale*democracy + authoritarian*democracy + socialdom*democracy,
                       data = johns_davies))
# --- table A1
screenreg(mod_full)

# substantively similar interaction between additive ethnocentrism and democracy treatment
summary(lm(SPEC_W1_E1q1 ~ democracy + Christian + SPEC_W1_E1casuals + socialdom + 
             Nationalism + authoritarian + gender + ConID + LibID + milpac + age65plus + 
             ethno_scale_additive*democracy + authoritarian*democracy + socialdom*democracy,
           data = johns_davies))

# as footnoted, the results are robust to dichotomization of the DV, i.e. the same DV representation used by Tomz and Weeks
johns_davies$SPEC_W1_E1q1_binary <- ifelse(johns_davies$SPEC_W1_E1q1 %in% c(4,5,6), 1, 0)
summary(glm(SPEC_W1_E1q1_binary ~ democracy + Christian + SPEC_W1_E1casuals + socialdom + 
             Nationalism + authoritarian + gender + ConID + LibID + milpac + age65plus + 
             ethno_scale*democracy + authoritarian*democracy + socialdom*democracy,
            family = "binomial",
           data = johns_davies))

# --- interaction plots
ethno_predict <- ggpredict(mod_full, c("ethno_scale", "democracy")) 
autho_predict <- ggpredict(mod_full, c("authoritarian", "democracy")) 
dom_predict <- ggpredict(mod_full, c("socialdom", "democracy")) 
plot_df <- 
  rbind(data.frame(strike = ethno_predict$predicted,
                   ethno = ethno_predict$x,
                   democ = ethno_predict$group,
                   ci_high = ethno_predict$conf.high,
                   ci_low = ethno_predict$conf.low,
                   model_iv = "ethnocentrism"),
        data.frame(strike = autho_predict$predicted,
                   ethno = autho_predict$x,
                   democ = autho_predict$group,
                   ci_high = autho_predict$conf.high,
                   ci_low = autho_predict$conf.low,
                   model_iv = "authoritarianism"),
        data.frame(strike = dom_predict$predicted,
                   ethno = dom_predict$x,
                   democ = dom_predict$group,
                   ci_high = dom_predict$conf.high,
                   ci_low = dom_predict$conf.low,
                   model_iv = "dominance"))

# --- figure A1
ggplot(subset(plot_df, model_iv == "ethnocentrism"), aes(x = ethno, y = strike, group = democ)) +
  geom_ribbon(aes(ymax = ci_high, ymin = ci_low), alpha = .1) +
  geom_line(aes(color = democ), linewidth = 2) +
  scale_color_manual(labels = c("No", "Yes"), values=c("#92978a","gray30",  "#c5bec2")) +
  labs(y = "Support for Strike", x = "Ethnocentrism", color = "Democracy?") +
  theme_minimal() +
  scale_y_continuous(n.breaks = 6) +
  theme(legend.title = element_text(size = 16),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 18, margin = margin(t=10)),
        axis.title.y = element_text(size = 18, margin = margin(r=10)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        legend.key = element_rect(color = NA, fill = NA),
        legend.key.size = unit(1, "cm"))

# --- figure A2
ggplot(subset(plot_df, model_iv == "dominance"), aes(x = ethno, y = strike, group = democ)) +
  geom_ribbon(aes(ymax = ci_high, ymin = ci_low), alpha = .1) +
  geom_line(aes(color = democ), linewidth = 2) +
  scale_color_manual(labels = c("No", "Yes"), values=c("#92978a","gray30",  "#c5bec2")) +
  labs(y = "Support for Strike", x = "Social Dominance Orientation", color = "Democracy?") +
  theme_minimal() +
  scale_y_continuous(n.breaks = 6) +
  coord_cartesian(ylim = c(0,3.5)) +
  theme(legend.title = element_text(size = 16),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 18, margin = margin(t=10)),
        axis.title.y = element_text(size = 18, margin = margin(r=10)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        legend.key = element_rect(color = NA, fill = NA),
        legend.key.size = unit(1, "cm"))

# --- figure A3
ggplot(subset(plot_df, model_iv == "authoritarianism"), aes(x = ethno, y = strike, group = democ)) +
  geom_ribbon(aes(ymax = ci_high, ymin = ci_low), alpha = .1) +
  geom_line(aes(color = democ), linewidth = 2) +
  scale_color_manual(labels = c("No", "Yes"), values=c("#92978a","gray30",  "#c5bec2")) +
  labs(y = "Support for Strike", x = "Authoritarianism", color = "Democracy?") +
  theme_minimal() +
  scale_y_continuous(n.breaks = 6) +
  coord_cartesian(ylim = c(0,3.5)) +
  theme(legend.title = element_text(size = 16),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 18, margin = margin(t=10)),
        axis.title.y = element_text(size = 18, margin = margin(r=10)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        legend.key = element_rect(color = NA, fill = NA),
        legend.key.size = unit(1, "cm"))

# --- conditional average treatment effects
johns_davies$ethnocentrism_binary <- ifelse(johns_davies$ethno_scale > median(na.omit(johns_davies$ethno_scale)), 1, 0)
# the democracy1 coef represents respondents at or below the ethno median;
# the democracy1:ethnocentrism_binary coef represents respondents above the ethno median
summary(mod_full <- lm(SPEC_W1_E1q1 ~ democracy + Christian + SPEC_W1_E1casuals  + socialdom + ethnocentrism_binary*democracy +
                         Nationalism + authoritarian + gender + ConID + LibID + milpac + age65plus + authoritarian*democracy +
                         socialdom*democracy,
                       data = johns_davies))

